        Z=incomevals;
        capital =0;
        
  
        phi=phivec(Qrho);
        kinit=0;
        Nvill=vss+1;
        
        
        NAiy= 350;              % number of grid points for kprime
        grid = 2;               % choose how you want to space the grid:
        NZ = size(Z,1);
        Nstar=NZ;
        nAiy = ones(1,Nstar); nnAiy = ones(NAiy,1); epsilon = 1e-5;
       kmin =  max(phi,phinat); kmax=150;
        
        if grid == 1;                     % grid when equally spaced
                kp = linspace(kmin,kmax,NAiy)';
                
            elseif grid == 2;                 % grid when logarithmically spaced
                kptemp = linspace(0,log(kmax+1-kmin),NAiy)';
                kp = exp(kptemp)-1+kmin;
                
            elseif grid == 3;                 % Chebyshec collocation nodes.
                YY     = -cos((2*[1:NAiy]'-1)*pi/(2*NAiy));
                kp = (YY+1)*(kmax-kmin)/2+kmin;
            end
            
            % ======================================================================
            % IV.b Start loop over decision rules
            % ======================================================================
              m0 = 0;                         % Initial value for cash on hand.
            mp = nnAiy*Z'+(1+r)*kp*nAiy;             % Cash on hand tomorrow.
            if capital==1
                kpp = (1-delta)*kp*nAiy;                % Initial guess for tomorrow's policy function (here zero net-investment).
            else
                kpp= kp*nAiy;
            end
            d22=1;
            it=1;
             while d22>1e-10
                cp=max(mp-max(kpp,phi),1e-10); % take a low value of c where rule would imply negative consumption: Indada
                % conditions apply
                EMUp = (cp).^(-rho)*Q1';
                m = ((beta*(1+r)*EMUp).^(-1/rho))+kp*nAiy;                  % Cash on hand today.
                d22 = max(max(abs(m0-m)./(1+abs(m))));                   % Update convergence measure.
                m0 = m;                                                  % Update initial value for cash on hand
                for i=1:length(Z)
                    kpp(:,i) = (interp1(m(:,i),kp,mp(:,i),'linear','extrap'));    % Update the policy function by interpolation.
                    kpp(:,i) = max(kpp(:,i),phi);
                end
                it=it+1
             end
                 kppphi=max(kpp,phi*ones(size(kpp)));  
       
       